This is a tutorial that shows how to make a map of SARS-CoV-2 RNA sewage sampling results in the Netherlands. We’ll start by importing the data. The data comes from this website (Description only provided in Dutch😢) : https://data.rivm.nl/meta/srv/chi/catalog.search#/metadata/a2960b68-9d3f-4dc3-9485-600570cd52b9?tab=relations

For a more general explanation of SARS-CoV2 sewage sampling in the Netherlands (In English!😊): https://www.rivm.nl/en/covid-19/sewage

Load the sewage RNA data

#make sure you've got the right working directory!
# getwd()    and   setwd()

nld_sewageRNA <- read.csv("nld.covid.sewage.csv")

#these are spatial packages we'll need to load
library(sf)
library(GISTools)
library(rgdal)
library(sp)
library(tmap)

#this is how you load the Netherlands shape file
# It's important that all the other .dbr, .prj files are stored in the same directory as the .shp file

nld.sf <- st_read("NLD_adm2.shp")
## Reading layer `NLD_adm2' from data source 
##   `C:\Users\rshea\Desktop\old computer\Thesis code\NLD_adm2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 491 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 3.360782 ymin: 50.75517 xmax: 7.29271 ymax: 53.55458
## Geodetic CRS:  WGS 84
#here's what the shape file looks like

plot(nld.sf$geometry)

Convert to RNA sewage sampling data to a spatial format

Next we will convert the coordinates from the excel sheet into a spatial object in the sp package and use the rest of the data from the spreadsheet to make a SpatialPointsData Frame

RNA.coords <- data.frame(x=nld_sewageRNA$X_coordinate, y=nld_sewageRNA$Y_coordinate)

nld_sewageRNA.sp <- SpatialPointsDataFrame(RNA.coords, nld_sewageRNA)

Reproject the SpatialPointsDataFrame

Now we need to assign a projection to the spatial points in order to map. The website that had the sewage sampling data provided an epsg code EPSG: 28992

Which you can then use to look up a corresponding proj4 code for on this website: https://spatialreference.org/ref/epsg/?search=28992&srtext=Search

There’s also this website that provides more info about the process:

https://www.earthdatascience.org/courses/earth-analytics/spatial-data-r/understand-epsg-wkt-and-other-crs-definition-file-types/

Then you can use the proj4string(object)<- “+[your proj4string]” command to set the correct projection

proj4string(nld_sewageRNA.sp)<- "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs "
## Warning in showSRID(uprojargs, format = "PROJ",
## multiline = "NO", prefer_proj = prefer_proj): Discarded
## datum Unknown based on Bessel 1841 ellipsoid in Proj4
## definition
#then we can convert the sp nld_sewageRNA.sp points to simple features (sf) 
#I'm not sure if this step is actually necessary...,but generally speaking sf is easier to work with....but unfortunately also newer so some earlier spatial methods haven't yet been adapted to work in sf, so you still have to run things in sp  

nld_sewageRNA.sf <- st_as_sf(nld_sewageRNA.sp)

Map the sewage sampling stations

For the sewage sampling stations mapping we use thematic mapping and we set the mapping mode to view instead of plot. That’s the step that makes an interactive webmap instead of just a static picture

This new map makes the following changes

  1. Dots are proportional in size to the amount of RNA in the sample

  2. The color of the dot corresponds to the date of the sample

  3. The bubbles were made semi-transparent using the alpha term

#this creates a separate layer for names which gives you more choices for things to turn and off in the viewer
nld.admin.labels <- nld.sf

#if this is changed to "plot" you get a static map
tmap_mode("view")
## tmap mode set to interactive viewing
#this creates the admin regions
tm_shape(nld.sf)+
  tm_polygons(col="purple", alpha = 0.05)+
  #alpha controls the transparency of the polygons
  
  #this adds the sampling result points
  tm_shape(nld_sewageRNA.sf)+
  tm_bubbles("RNA_per_ml",col="Date_measurement", alpha=0.5)+
  
  #this creates the names layer
  tm_shape(nld.admin.labels)+
  tm_text("NAME_2", size=1)
## Warning: Number of levels of the variable
## "Date_measurement" is 170, which is larger than
## max.categories (which is 30), so levels are combined.
## Set tmap_options(max.categories = 170) in the layer
## function to show all levels.
## Legend for symbol sizes not available in view mode.
#if this is changed to "plot" you get a static map
tmap_mode("view")
## tmap mode set to interactive viewing
#this creates the admin regions
sewageRNA.map <- tm_shape(nld.sf)+
  tm_polygons(col="purple", alpha = 0.05)+
  #alpha controls the transparency of the polygons
  
  #this adds the sewage sampling points
  tm_shape(nld_sewageRNA.sf)+
  tm_bubbles("RNA_per_ml",
             col="Date_measurement", 
             palette= "PuBuGn",
             alpha=0.5)+
  
  #this creates the names layer
  tm_shape(nld.admin.labels)+
  tm_text("NAME_2", size=1)

sewageRNA.map
## Warning: Number of levels of the variable
## "Date_measurement" is 170, which is larger than
## max.categories (which is 30), so levels are combined.
## Set tmap_options(max.categories = 170) in the layer
## function to show all levels.
## Legend for symbol sizes not available in view mode.

Here are some useful documentation files:

https://www.rdocumentation.org/packages/tmap/versions/3.0/topics/tmap_animation

https://www.rdocumentation.org/packages/tmap/versions/3.0/topics/tm_facets

This step merges all the boundaries in the shape file together. We’ll need this file for making our GIF

netherlands.outline.sf <- st_union(nld.sf)
# the code is set to not run to save space in the tutorial, change eval to TRUE or copy and paste this code into R to make the gif

tmap_mode("plot")

netherlands.outline.sf <- st_union(nld.sf)

Nld_sewageRNA.gif <-  tm_shape(netherlands.outline.sf)+
  tm_polygons(col="purple", alpha = 0.05)+
  tm_layout(title = "SARS-CoV2 RNA in sewage at various sites in the Netherlands")+
  #alpha controls the transparency of the polygons
  
  #this adds the sewage sampling points
  tm_shape(nld_sewageRNA.sf)+
  tm_bubbles("RNA_per_ml",
             col="red",
             alpha=0.8)+
  tm_facets(along = "Date_measurement" , sync = FALSE, as.layers = TRUE)

library(gifski)

tmap_animation(Nld_sewageRNA.gif, filename="Nld_sewageRNA.gif", delay = 100, restart.delay = 300 )

Finally save the file as “NLD COVID19 sewage sampling.Rmd” in your working directory. Then load the markdown library using library(rmarkdown) then use the command render(“NLD COVID19 sewage sampling.Rmd”) to create the HTML page in your library

This code can be used to save a full screen version of the map:

?tmap_save

tmap_save(sewageRNA.map)